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A study of QCD at non-zero chemical potential, p, and temperature, T, is per- 
formed using the lattice technique. The transition temperature (between the con- 
fined and dcconfined phases) is determined as a function of pi and is found to be in 
agreement with other work. In addition the variation of the pressure and energy 
density with p is obtained for p > 0. These results are of particular relevance for 
heavy-ion collision experiments. 



1. Introduction 

The QCD phase diagram has come under increasing experimental and theo- 
retical scrutiny over the last few years. On the experimental side, very recent 
studies of compact astronomical objects have suggested that their cores con- 
tain "quark matter" , i.e. QCD in a new, unconfmcd phase where the basic 
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units of matter are quarks, rather than nuclei or nuclcons 1 . More terrestrially, 
heavy ion collision experiments, such as those performed at RHIC and CERN, 
are also believed to be probing unconfincd QCD 2 . On the theoretical side, the 
study of QCD under these extreme densities and temperatures has proceeded 
along several fronts. One of the most promising areas of research is the use of 
lattice techniques to study either QCD itself, or model theories which mimic 
the strong interaction 3 . Clearly the most satisfying approach would be the 
former, i.e. a direct lattice study of QCD at various coordinates (T, fi) in its 
phase space (/x is the chemical potential for the quark number). However, until 
recently, this has proved intractable at a practical level for very fundamental 
reasons. This is because the Monte Carlo integration technique, which is at 
the heart of the (Euclidean) lattice approach, breaks down when zz ^ 0. This 
work summarises one new approach which overcomes this problem and has 
made progress for /j, ^ and T ^ 0. 

In the next section a summary is given of the lattice technique and the 
problem incurred when /i ^ 0. Section 3 describes the method used to overcome 
these difficulties and section 4 outlines the simulation details. The next two 
sections apply the method to variations in m and /i, and section 7 describes 
calculations of the pressure and energy density as functions of fi. 

A full account of this work is published elsewhere 4 . 



2. Lattice technique 

On the lattice, the quark fields, ip(x), are defined on the sites, x, and the 
gluonic fields, U^x), on the links x — > x + fx. Observables are then calculated 
via a Monte Carlo integration approach: 

<n> = — - — V' ny>,$,u) (1) 

where J2' represents a sum over configurations {Ujip,^} which are se- 
lected with probability proportional to the Boltzmann weight P({U, ip, V>}) oc 
e -s({u,ip,tp}) w h ere 5* i s a suitably defined (Euclidean) gauge-invariant action. 
The fermionic part of this action is 



Sjr = (#+m)#r). (2) 

= M 

For = it can be shown that this action produces a (real-valued) positive 
Boltzmann weight. 

Calculations at non-zero temperature, T ^ 0, can be performed by using 
a lattice with a finite temporal extent of N t a = 1/T, where N t is the number 
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of lattice sites in the time dimension. In practice, T is varied by changing the 
gauge coupling, go, and hence (through dimensional transmutation) the lattice 
spacing, a, rather than by changing N t (which can only be changed in discrete 
steps!). 

The chemical potential is introduced into the system via an additional term 
in the quark matrix M, proportional to the Dirac matrix, 70, 

M->M + ^ . (3) 

For fj, ^ 0, this leads to a complex-valued Boltzmann "weight" which can 
therefore no longer be used as a probability distribution, and, hence, the Monte 
Carlo integration procedure is no longer applicable. This is known as the Sign 
Problem and has plagued more than a decade of lattice calculations of QCD 
at fj, ^ 0. 



3. Reweighting 

This section outlines the Ferrenberg-Swcndsen reweighting approach 5 which is 
used to overcome the sign problem detailed in the previous section. Observ- 
ables at one set of parameter values (/3, m, (where (3 = 6/g%, and m is the 
quark mass) can be calculated using an ensemble generated at another set of 
parameters (/3o, m-o> Mo) as follows, 

< fi W.") = Z((i] m ^) j DUndetM(m,ric- s ^ (4) 

(Q, e ( lndet M ( m ^)- lndot M ( m o,Mo)) e -5 9 (/3)+S 9 (/3 )\ 

/ c (lndotM(m, AI )-lndotM(m ,Mo)) e -S's(/3)+S g (/3o)\ ' ^ ' 

\ / (0o,mo,no) 

Here S g is the gauge action. 

In principle, Eq.(5) can be used to map out the entire phase diagram of 
QCD. However, it has been found that its naive application fails since the 
relative size of fluctuations in both numerator and denominator tend to grow 
exponentially with the volume of the system studied. This is a signal that the 
overlap between the ensemble at (/3o, mo, Ho) is small compared with that at 
(/3,m,n). 

One study which has had success in using Eq.(5) is that of Fodor and Katz 6 . 
They apply the reweighting approach not at an arbitrary point, (T, fi), in the 
phase diagram, but rather they trace out the phase transition line T c (fi). This 
method is successful presumably because the overlap between the ensembles 
remains high along the "coexistence" line which defines the transition. 

This paper utilises an alternative approach which Taylor expands Eq.(5) as 
a function of \x (or m) and hence estimates the derivatives of various quantities 
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w.r.t. n (or m). Derivatives up to the second order are considered, thus for 
< 51 > (in the case of a Taylor expansion in /i) we have 

, ON {(go + "lM + ^2M 2 ) exp(i? lM + i? 2M 2 - AS 9 )) (/Wo) 

{ ){M ~ <ex P (i?^ + AS g )) {M ' W 

where fi = 0. In Eq.(6), R n is the n— th derivative of the fermionic reweighting 
factor in Eq.(5), i.e. 

/ dctM( M ) \ ^ M " 0"lndetM(O) _^ 
ln ^detM(0)J-2-^ = 2^" • (7) 

v / / n=l n=l 

The Q n are similarly the n— th derivatives of the observable fi. 

Two observables are studied: the chiral condensate < ipip > and the 
Polyakov Loop, L. Since L is a pure gluonic quantity, defined as 



1 Wt 

<i>=< FE Tr Il c/ *( f ^) > ' ( 8 ) 

x t— U 

all of its derivatives are zero. (Here, 14 is the spatial volume.) However, the 
expansion of < ipip > is more challenging since it is defined as 

(iPiP) ~ (trM- 1 ), (9) 

and hence the application of Eq.(6) requires determinations of d n M~ 1 /dfi n . 

The susceptibilities, \i °f both < ipip > and L are defined as usual by their 
fluctuations, e.g. 

XL = (volume factor) x (< L 2 > - < L > 2 ). (10) 

These susceptibilities have a maximum at the transition point, /3 C , and hence 
can be used to determine the transition point /3 c (m, /u). 



4. Simulation details 

The lattice calculations were performed using a "p4-improved" discretisation 
of the continuum action which is a sophisticated lattice action maintaining 
rotational invariance of the free fermion propagator up to 0(p 4 ) 7,8 (p is the 
momentum here). Two dynamical flavours of quarks were used with a 16 3 x 4 
lattice. Simulations were performed at quark mass, m = 0.1 and 0.2 which 
correspond to (unphysically heavy) pseudoscalar-vector meson mass ratios of 
Mps/My fts 0.70 and 0.85. Approximately 400,000 configurations were gener- 
ated in total using around 6 months of a 128-node (64Gflop peak) APEmillc 
in the University of Wales Swansea. 
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5. Results for mass reweighting 

As a check of our method, we first use mass reweighting, since, unlike the \i ^ 
case, there is no theoretical difficulty in simulating at virtually any value of to, 
and hence there are published data at a number of different to values readily 
available for comparison. Reweighting in quark mass is simply a matter of 
setting /x = /io = in Eq.(5) and Taylor expanding in to rather than n in 
Eqs.(6 & 7). (Note that dM/dm = 1 and d n M/dm n = for n > 1.) 

We use the peak position in both the chiral condensate and Polyakov Loop 
susceptibilities, x^ l *° determine the phase transition point (3 c (m). Figure 1 
illustrates this by plotting \^ as a function of (3 (for to = 0.2). This shows the 
variation in the peak position of as m is changed in steps of 0.01 around 
too = 0.2. 
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Figure 1. The quark mass dependency of Xtfj, as a function of f) in the neighbourhood of 
m = 0.2 for the Chiral Susceptibility. 



Once these peak positions have been determined (for both and Xl), 
(3 C can be plotted as a function of to and a comparison can be made with other 
determinations. This is done in Figure 2 where the f3 c values from earlier 
work 7 are also shown. The line segments around our data points at m = 
0.1 and 0.2 indicate the gradient of (3 C as a function of to using our Taylor- 
expanded reweighting technique. Shown in Figure 2 are results from both the 
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chiral condensate and Polyakov loop. These are both in perfect agreement (as 
expected). Furthermore, our method agrees with previously published data 
confirming the validity of the approach. 
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Figure 2. The transition "temperature", /3 C , as a function of m in comparison with previous 
work 7 



6. Results for fi reweighting 

We now turn to reweighting in chemical potential, /U. As in 6 , rather than 
applying the method to arbitrary parameter values, we trace out the transition 
point P c {ix). Using Eq.(6), < ipip > and the Polyakov Loop, L, are calculated 
as a function of n together with their susceptibilities, \^ L (see Eq.(10)). In 
Figure 3, we plot against [3 for various (i. Note that the peak position 
moves as [i changes. The determination of the transition point (3 C from both 
the chiral condensate and the Polyakov loop (not shown here) are found to be 
in agreement. 

Because we have calculated all quantities to 0(/j, 2 ) we can extract f} c and 
fit it to a quadratic in [i. (In fact, it can be shown 4 that the first derivative 
d/? c /d/z = 0.) We find d 2 /3 c /d/i 2 = 1.1 ± 50% for both quark masses m = 0.1 
and 0.2. 
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to convert d 2 /3 c /d/x 2 into physical units, with the beta-function d/3/da coming 
from string tension data 7 . We find T c ^f « —0.14 at ma = 0.1. 

Figure 4 shows the phase transition curve T c (/x) obtained from this method. 
On this graph we have plotted the Fodor-Katz point 6 which is within our errors, 
confirming our method. Also shown is the /i value corresponding to RHIC. It 
is interesting to extrapolate the curve T c (//) to the horizontal axis (as shown). 
It is known that the transition (at T«0) between ordinary hadronic matter 
and quark matter occurs at around /x w 400 MeV. This is at a smaller value 
of \x than the horizontal intercept of our data indicating (not surprisingly) the 
presence of higher order terms in the Taylor expansion and/or a breakdown in 
our method at these large values of fx. 

This motivates the question: for what range in /x do we expect our method 
to be accurate (and converge to the correct answer)? We have studied this 
issue by calculating the complex phase, 9, of the fermionic determinant (which 
enters in the reweighting factor in Eq.(5)), i.e. 

detM= |detM|e ie (12) 
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Figure 4. The transition temperature, T c , as a function of fi. The diamond symbol is the 
endpoint of the first order transition obtained by Fodor and Katz 6 

The reweighting method will fail when the fluctuations (standard deviation) 
in 9 are larger than 0(n/2). Taylor expanding Eq.(12) and noting that only 
odd derivatives contribute to the complex phase 9, we find that the standard 
deviation A9 ~ 0(n/2) at around n/T c ~ 0.5. Since this is around five times 
the [i value of RHIC, we can confirm that our method is applicable for RHIC 
physics. 

An interesting dynamical quark effect can be uncovered when studying the 
Polyakov Loop susceptibility, \l- For /i < 0, anti-quarks are dynamically 
generated which screen colour charge. This leads to a reduction in the free 
energy of a single quark, and a corresponding reduction in the strength of the 
singularity at the transition. We observe this effect by noting that the peak 
height of xl is smaller for /i < compared with /i > 4 . 

7. Pressure and energy density 

Of great interest for heavy ion collision experiments is the study of the pressure, 
p, and energy density e and their /i dependence. We can obtain estimates of p 
by employing the integral method : 



P 




(13) 
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The first derivative of p w.r.t. ji is related to the quark number density, 

TdlnZ 

n «=v^ (14) 

and the second derivative to the singlet quark number susceptibility, \S- Both 
n q and xs can be calculated in terms of the quark matrix, M. 

Using the above to estimate p at the RHIC point we find that p increases 
by around 1% from its /j, = value. 

The energy density, e, can be obtained from 



e — 3p 1 d In Z f 1 da 



T 4 VT 3 dp \adp 



(15) 



The derivatives of Eq.(15) can be expressed in terms of n q and \s- Combining 
this with the above calculation of p, we obtain a value for e alone. We find 
that, at the RHIC point, there is again only a 1% deviation from e(yu = 0). 

Finally we study the variation of p and e along the transition line T c (/x). 
(The above calculations were performed at fixed T.) Our aim is to determine 
whether p and e are constant along T c (n). The constant p line is defined as 

with a similar expression for the constant e line. Using the above and the value 
determined earlier for the rate of change of T c with n, we find that the value 
of both p and e along the transition line T c (n) is consistent with zero with our 
current precision 4 . 



8. Conclusions 

This work (which is published in full elsewhere 4 ) has outlined a new method of 
determining thermal properties of QCD at non-zero chemical potential from the 
lattice. This approach is based on Taylor expanding the Ferrenberg-Swendsen 
reweighting scheme. 

Using this method, the susceptibilities in the chiral condensate and 
Polyakov loop were determined and their peak positions used to define the 
transition point, T c . As a warmup exercise, the reweighting technique was 
used to determine the transition point as a function of the quark mass, m 
confirming earlier work. The method was then applied to obtaining T c as a 
function of chemical potential, ^i, confirming the work of Fodor and Katz 6 . 
Very recent work of de Forcrand and Philipsen, who studied the transition 
temperature for imaginary fi and then analytically continued these results to 
real- valued [i also confirm our results 10 . 
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The region of applicability of our method was studied by calculating the 
fluctuations in the phase of the reweighting factor. This region was found to be 
substantial and easily covers the physically interesting values of fi appropriate 
for RHIC physics. 

We also extracted information about the pressure, p, and energy density, 
e, as a function of chemical potential. We found that the variation in these 
quantities from their values at ji — is tiny. This leads us to conclude that 
RHIC physics well approximated by [i = physics. Furthermore we find that 
p and e are approximately constant along the transition line T c (fi). 

The success of this work motivates the use of lighter, more physical quark 
masses, and the study of (2+1) dynamical flavours to correctly model real 
world physics. 

Acknowledgments 

The authors would like to acknowledge the Particle Physics and Astronomy 
Research Council for the award of the grant PPA/G/S/1999/00026 and the 
European Union for the grant ERBFMRX-CT97-0122. CRA would like to 
thank the Department of Mathematics, University of Queensland, Australia 
for their kind hospitality while part of this work was performed. 

References 

1. see http : //wwwl . msf c . nasa . gov/NEWSR00M/news/releases/2002/02-082 . html 

(NASA Press release); and J. Drake et al., The Astrophysical Journal, June 20 
2002 to appear. 

2. For a recent announcement from CERN, see 
http://cern.web. cern. ch/CERN/Announcements/2000/NewStateMatter/. 

3. For a review, see S.J. Hands, Nucl.Phys. Proc.Suppl. 106 (2002) 142, 
hep-lat/0109034. 

4. C.R. Allton, S. Ejiri, S.J. Hands, O. Kaczmarek, F. Karsch, E. Laermann, Ch. 
Schmidt, L. Scorzato, hep-lat/0204010. 

5. A.M. Ferrenberg and R.H. Swendsen, Phys. Rev. Lett. 61 (1988) 2635; Phys. Rev. 
Lett. 63 (1989) 1195. 

6. Z. Fodor and SD. Katz, hep-lat/0104001; JHEP 0203:014,2002 
hep-lat/0106002. 

7. F. Karsch, E. Laermann, and A. Peikert, Phys. Lett. B478 (2000) 447; Nucl. Phys. 
B605 (2001) 579. 

8. U.M. Heller, F. Karsch, and B. Sturm, Phys. Rev. D60 (1999) 114502. 

9. J. Engels, J. Fingberg, F. Karsch, D. Miller and M. Weber, Phys. Lett. B252 
(1990) 625. 

10. Ph. de Forcrand and O. Philipsen, hep-lat/0205016. 



